1 Setup

suppressPackageStartupMessages(library(data.table))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(gplots))
suppressPackageStartupMessages(library(here))
suppressPackageStartupMessages(library(hyperSpec))
suppressPackageStartupMessages(library(RColorBrewer))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(VennDiagram))
suppressMessages(source(here("UPSCb-common/src/R/featureSelection.R")))
suppressMessages(source(here("UPSCb-common/src/R/plotMA.R")))
suppressMessages(source(here("UPSCb-common/src/R/volcanoPlot.R")))
suppressMessages(source(here("UPSCb-common/src/R/gopher.R")))
pal=brewer.pal(8,"Dark2")
hpal <- colorRampPalette(c("blue","white","red"))(100)
mar <- par("mar")

Load the genes of interests (goi)

gois <- list(
    myb.goi=scan(here("doc/MYB-goi.txt"),what="character"),
    mybr.goi=scan(here("doc/MYB-related-goi.txt"),what="character"),
    wrky.goi=scan(here("doc/WRKY-goi.txt"),what="character"),
    auxin.goi=scan(here("doc/auxin-related-goi.txt"),what="character"),
    pectin.goi=scan(here("doc/pectin-related-goi.txt"),what="character"))
  1. plot specific gene expression
"line_plot" <- function(dds=dds,vst=vst,gene_id){
    sel <- grepl(gene_id,rownames(vst))
    stopifnot(sum(sel)==1)

    
    p <- ggplot(bind_cols(as.data.frame(colData(dds)),
                                 data.frame(value=vst[sel,])),
               aes(x=Time,y=value,col=Experiment,group=Experiment)) +
            geom_point() + geom_smooth() +
            scale_y_continuous(name="VST expression") + 
            ggtitle(label=paste("Expression for: ",gene_id))
    suppressMessages(suppressWarnings(plot(p)))
    return(NULL)
}
  1. extract the DE results. Default cutoffs are from Schurch et al., RNA, 2016
"extract_results" <- function(dds,vst,contrast,
                              padj=0.01,lfc=0.5,
                              plot=TRUE,verbose=TRUE,
                              export=TRUE,default_dir="analysis/DE",
                              default_prefix="DE-",
                              labels=colnames(dds),
                              sample_sel=1:ncol(dds)){
    
    if(length(contrast)==1){
        res <- results(dds,name=contrast)
    } else {
        res <- results(dds,contrast=contrast)
    }
    
    if(plot){
        par(mar=c(5,5,5,5))
        volcanoPlot(res)
        par(mar=mar)
    }
    
    sel <- res$padj <= padj & abs(res$log2FoldChange) >= lfc & ! is.na(res$padj)
    
    if(verbose){
        message(sprintf("There are %s genes that are DE",sum(sel)))
    }
            
    if(export){
        if(!dir.exists(default_dir)){
            dir.create(default_dir,showWarnings=FALSE,recursive=TRUE,mode="0771")
        }
        write.csv(res,file=file.path(default_dir,paste0(default_prefix,"results.csv")))
        write.csv(res[sel,],file.path(default_dir,paste0(default_prefix,"genes.csv")))
    }
    if(plot){
        heatmap.2(t(scale(t(vst[sel,sample_sel]))),
                  distfun = pearson.dist,
                  hclustfun = function(X){hclust(X,method="ward.D2")},
                  trace="none",col=hpal,labRow = FALSE,
                  labCol=labels[sample_sel]
        )
    }
    return(list(all=rownames(res[sel,]),
                up=rownames(res[sel & res$log2FoldChange > 0,]),
                dn=rownames(res[sel & res$log2FoldChange < 0,])))
}

2 Laccaria bicolor

load(here("data/analysis/salmon/Lacbi-all-dds.rda"))

2.1 Normalisation for visualisation

vsd <- varianceStabilizingTransformation(dds,blind=FALSE)
vst <- assay(vsd)
vst <- vst - min(vst)
  • Genes from the GOI
rownames(vst) <- substr(rownames(vst),12,17)
  • Genes from the GOI
lapply(gois, function(goi){
        pdf(paste0(goi,"-Lacbi-candidates-line-plot.pdf"),width=12,height=8)
        par(mfrow=c(2,3))
        lapply(goi,function(x){
        if ((x %in% rownames(vst))){
            line_plot(dds,vst,x)
        }
    })
    dev.off()
})
## $myb.goi
## png 
##   2 
## 
## $mybr.goi
## png 
##   2 
## 
## $wrky.goi
## png 
##   2 
## 
## $auxin.goi
## png 
##   2 
## 
## $pectin.goi
## png 
##   2

2.2 Differential Expression

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
  • Dispersion estimation The dispersion estimation is adequate
plotDispEsts(dds)

The model used is:

Experiment * Time meaning that the Experiment and Time variable as well as their interaction Experiment:Time is considered. Because we cannot assume that these two variables explain all the variance in the data, there is also an Intercept for the linear model. This also implies that the model assumes ECM at 3 hours to be the baseline; i.e. everything is compared against it.

resultsNames(dds)
##  [1] "Intercept"             "Experiment_FLM_vs_ECM" "Time_7_vs_3"          
##  [4] "Time_14_vs_3"          "Time_21_vs_3"          "Time_28_vs_3"         
##  [7] "ExperimentFLM.Time7"   "ExperimentFLM.Time14"  "ExperimentFLM.Time21" 
## [10] "ExperimentFLM.Time28"

2.3 Results

In the following we look at the interaction specific genes; i.e. genes that changes at a given time transition in between experiments ### FLM vs. ECM at T3

Lb_3 <- extract_results(dds,vst,"Experiment_FLM_vs_ECM",
                        default_prefix="Labic_FLM-vs-ECM_T3_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==3)
## Loading required package: LSD
## There are 2100 genes that are DE

2.3.1 FLM vs. ECM at T7

Here we want to conmbine the effect of FLM-ECM at time T3 and the specific FLM:T7 interaction

Lb_7 <- extract_results(dds,vst,c(0,1,0,0,0,0,1,0,0,0),
                       default_prefix="Labic_FLM-vs-ECM_T7_",
                       labels=paste0(colData(dds)$Experiment,
                                     colData(dds)$Time),
                       sample_sel=colData(dds)$Time==7)
## There are 1222 genes that are DE

2.3.2 FLM vs. ECM at T14

Lb_14 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,1,0,0),
                   default_prefix="Labic_FLM-vs-ECM_T14_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==14)
## There are 362 genes that are DE

2.3.3 FLM vs. ECM at T21

Lb_21 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,1,0),
                        default_prefix="Labic_FLM-vs-ECM_T21_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==21)
## There are 381 genes that are DE

2.3.4 FLM vs. ECM at T28

Lb_28 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,0,1),
                        default_prefix="Labic_FLM-vs-ECM_T28_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==28)
## There are 673 genes that are DE

2.3.5 Venn Diagram

2.3.5.1 All DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Lb_3$all,
                            T7=Lb_7$all,
                            T14=Lb_14$all,
                            T21=Lb_21$all,
                            T28=Lb_28$all),
                       NULL,
                       fill=pal[1:5]))

2.3.5.2 UP DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Lb_3$up,
                            T7=Lb_7$up,
                            T14=Lb_14$up,
                            T21=Lb_21$up,
                            T28=Lb_28$up),
                       NULL,
                       fill=pal[1:5]))

2.3.5.3 DOWN DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Lb_3$dn,
                            T7=Lb_7$dn,
                            T14=Lb_14$dn,
                            T21=Lb_21$dn,
                            T28=Lb_28$dn),
                       NULL,
                       fill=pal[1:5]))

res.list <- list(Lb_3=list(all=substr(Lb_3$all,12,17),
                           up=substr(Lb_3$up,12,17),
                           dn=substr(Lb_3$dn,12,17)),
                 Lb_7=list(all=substr(Lb_7$all,12,17),
                           up=substr(Lb_7$up,12,17),
                           dn=substr(Lb_7$dn,12,17)),
                 Lb_14=list(all=substr(Lb_14$all,12,17),
                            up=substr(Lb_14$up,12,17),
                            dn=substr(Lb_14$dn,12,17)),
                 Lb_21=list(all=substr(Lb_21$all,12,17),
                            up=substr(Lb_21$up,12,17),
                            dn=substr(Lb_21$dn,12,17)),
                 Lb_28=list(all=substr(Lb_28$all,12,17),
                            up=substr(Lb_28$up,12,17),
                            dn=substr(Lb_28$dn,12,17)))

2.3.6 Gene Ontology enrichment

background <- rownames(vst)[featureSelect(vst,dds$Experiment,exp=0.1)]

enr.list <- lapply(res.list,function(r){
    lapply(r,gopher,background=background,task="go",url="lacbi2")
})
## Loading required package: jsonlite
## 
## Attaching package: 'jsonlite'
## The following object is masked from 'package:purrr':
## 
##     flatten
## No enrichments found in task: go
## No enrichments found in task: go
## No enrichments found in task: go
## No enrichments found in task: go
dev.null <- lapply(names(enr.list),function(n){
    r <- enr.list[[n]]
    if (! is.null(r$all$go)){
        write_delim(r$all$go,path=file.path(file.path(here("data/analysis/DE",
                                                       paste0(n,"-all-DE-genes_GO-enrichment.txt")))))
        write_delim(r$all$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                        paste0(n,"-all-DE-genes_GO-enrichment_for-REVIGO.txt")))))
    }
    if (! is.null(r$up$go)){
        write_csv(r$up$go,path=file.path(file.path(here("data/analysis/DE",
                                                    paste0(n,"-up-DE-genes_GO-enrichment.txt")))))
        write_delim(r$up$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                       paste0(n,"-up-DE-genes_GO-enrichment_for-REVIGO.txt")))))    
    }
    if (! is.null(r$dn$go)){
        write_csv(r$dn$go,path=file.path(file.path(here("data/analysis/DE",
                                                    paste0(n,"-down-DE-genes_GO-enrichment.txt")))))
        write_delim(r$dn$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                       paste0(n,"-down-DE-genes_GO-enrichment_for-REVIGO.txt")))))    
    }
})

3 Populus tremula

load(here("data/analysis/salmon/Potra-104-127-139-removed-dds.rda"))

3.1 Normalisation for visualisation

vsd <- varianceStabilizingTransformation(dds,blind=FALSE)
vst <- assay(vsd)
vst <- vst - min(vst)

3.2 Gene of interest

  • Potra000962g07909 The gene is not expressed
line_plot(dds,vst,"Potra000962g07909")

## NULL
  • Genes from the GOI
lapply(gois, function(goi){
    pdf(paste0(goi,"-Potra-candidates-line-plot.pdf"),width=12,height=8)
    par(mfrow=c(2,3))
    lapply(goi,function(x){
        if ((x %in% rownames(vst))){
            line_plot(dds,vst,x)
        }
    })
    dev.off()
})
## $myb.goi
## png 
##   2 
## 
## $mybr.goi
## png 
##   2 
## 
## $wrky.goi
## png 
##   2 
## 
## $auxin.goi
## png 
##   2 
## 
## $pectin.goi
## png 
##   2

3.3 Differential Expression

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
  • Dispersion estimation The dispersion estimation is adequate
plotDispEsts(dds)

The model used is:

Experiment * Time meaning that the Experiment and Time variable as well as their interaction Experiment:Time is considered. Because we cannot assume that these two variables explain all the variance in the data, there is also an Intercept for the linear model. This also implies that the model assumes Cont at 3 hours to be the baseline; i.e. everything is compared against it.

resultsNames(dds)
##  [1] "Intercept"              "Experiment_ECM_vs_Cont" "Time_7_vs_3"           
##  [4] "Time_14_vs_3"           "Time_21_vs_3"           "Time_28_vs_3"          
##  [7] "ExperimentECM.Time7"    "ExperimentECM.Time14"   "ExperimentECM.Time21"  
## [10] "ExperimentECM.Time28"

3.4 Results

In the following we look at the interaction specific genes; i.e. genes that changes at a given time transition in between experiments ### ECM vs. Cont at T3

Pa_3 <- extract_results(dds,vst,"Experiment_ECM_vs_Cont",
                        default_prefix="Potra_ECM-vs-Cont_T3_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==3)
## There are 230 genes that are DE

3.4.1 ECM vs. Cont at T7

Pa_7 <- extract_results(dds,vst,c(0,1,0,0,0,0,1,0,0,0),
                        default_prefix="Potra_ECM-vs-Cont_T7_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==7)
## There are 205 genes that are DE

3.4.2 ECM vs. Cont at T14

Pa_14 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,1,0,0),
                         default_prefix="Potra_ECM-vs-Cont_T14_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==14)
## There are 76 genes that are DE

3.4.3 ECM vs. Cont at T21

Pa_21 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,1,0),
                         default_prefix="Potra_ECM-vs-Cont_T21_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==21)
## There are 1276 genes that are DE

3.4.4 ECM vs. Cont at T28

Pa_28 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,0,1),
                         default_prefix="Potra_ECM-vs-Cont_T28_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==28)
## There are 416 genes that are DE

3.4.5 Venn Diagram

3.4.5.1 All DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$all,
                            T7=Pa_7$all,
                            T14=Pa_14$all,
                            T21=Pa_21$all,
                            T28=Pa_28$all),
                       NULL,
                       fill=pal[1:5]))

3.4.5.2 UP DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$up,
                            T7=Pa_7$up,
                            T14=Pa_14$up,
                            T21=Pa_21$up,
                            T28=Pa_28$up),
                       NULL,
                       fill=pal[1:5]))

3.4.5.3 DOWN DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$dn,
                            T7=Pa_7$dn,
                            T14=Pa_14$dn,
                            T21=Pa_21$dn,
                            T28=Pa_28$dn),
                       NULL,
                       fill=pal[1:5]))

res.list <- list(Pa_3=Pa_3,
                 Pa_7=Pa_7,
                 Pa_14=Pa_14,
                 Pa_21=Pa_21,
                 Pa_28=Pa_28)

3.4.6 Gene Ontology enrichment

background <- rownames(vst)[featureSelect(vst,dds$Experiment,exp=0.1)]

enr.list <- lapply(res.list,function(r){
    lapply(r,gopher,background=background,task="go",url="potra")
})
## No enrichments found in task: go
## No enrichments found in task: go
## No enrichments found in task: go
dev.null <- lapply(names(enr.list),function(n){
    r <- enr.list[[n]]
    if (! is.null(r$all$go)){
        write_delim(r$all$go,path=file.path(file.path(here("data/analysis/DE",
                                                           paste0(n,"-all-DE-genes_GO-enrichment.txt")))))
        write_delim(r$all$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                            paste0(n,"-all-DE-genes_GO-enrichment_for-REVIGO.txt")))))
    }
    if (! is.null(r$up$go)){
        write_csv(r$up$go,path=file.path(file.path(here("data/analysis/DE",
                                                        paste0(n,"-up-DE-genes_GO-enrichment.txt")))))
        write_delim(r$up$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                           paste0(n,"-up-DE-genes_GO-enrichment_for-REVIGO.txt")))))    
    }
    if (! is.null(r$dn$go)){
        write_csv(r$dn$go,path=file.path(file.path(here("data/analysis/DE",
                                                        paste0(n,"-down-DE-genes_GO-enrichment.txt")))))
        write_delim(r$dn$go[,c("id","padj")],path=file.path(file.path(here("data/analysis/DE",
                                                                           paste0(n,"-down-DE-genes_GO-enrichment_for-REVIGO.txt")))))    
    }
})

4 Populus trichocarpa

load(here("data/analysis/salmon/Potri-all-dds.rda"))

4.1 Normalisation for visualisation

vsd <- varianceStabilizingTransformation(dds,blind=FALSE)
vst <- assay(vsd)
vst <- vst - min(vst)
  • Potri.006G134800 The gene is not expressed
line_plot(dds,vst,"Potri.006G134800")

## NULL
  • Potri.006G221800 The gene is medium expressed and show time-lagged expression patterns between the experiments
line_plot(dds,vst,"Potri.006G221800")

## NULL
  • Potri.002G173900 The gene is medium expressed and show initially opposing expression patterns between experiments
line_plot(dds,vst,"Potri.002G173900")

## NULL
  • Potri.004G088100 The gene his very lowly expressed and shows somewhat opposite expression patterns
line_plot(dds,vst,"Potri.004G088100")

## NULL

4.2 Differential Expression

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
  • Dispersion estimation The dispersion estimation is adequate
plotDispEsts(dds)

The model used is:

Experiment * Time meaning that the Experiment and Time variable as well as their interaction Experiment:Time is considered. Because we cannot assume that these two variables explain all the variance in the data, there is also an Intercept for the linear model. This also implies that the model assumes Cont at 3 hours to be the baseline; i.e. everything is compared against it.

resultsNames(dds)
##  [1] "Intercept"              "Experiment_ECM_vs_Cont" "Time_7_vs_3"           
##  [4] "Time_14_vs_3"           "Time_21_vs_3"           "Time_28_vs_3"          
##  [7] "ExperimentECM.Time7"    "ExperimentECM.Time14"   "ExperimentECM.Time21"  
## [10] "ExperimentECM.Time28"

4.3 Results

In the following we look at the interaction specific genes; i.e. genes that changes at a given time transition in between experiments ### ECM vs. Cont at T3

Pa_3 <- extract_results(dds,vst,"Experiment_ECM_vs_Cont",
                        default_prefix="Potri_ECM-vs-Cont_T3_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==3)
## There are 311 genes that are DE

4.3.1 ECM vs. Cont at T7

Pa_7 <- extract_results(dds,vst,c(0,1,0,0,0,0,1,0,0,0),
                        default_prefix="Potri_ECM-vs-Cont_T7_",
                        labels=paste0(colData(dds)$Experiment,
                                      colData(dds)$Time),
                        sample_sel=colData(dds)$Time==7)
## There are 219 genes that are DE

4.3.2 ECM vs. Cont at T14

Pa_14 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,1,0,0),
                         default_prefix="Potri_ECM-vs-Cont_T14_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==14)
## There are 220 genes that are DE

4.3.3 ECM vs. Cont at T21

Pa_21 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,1,0),
                         default_prefix="Potri_ECM-vs-Cont_T21_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==21)
## There are 1899 genes that are DE

4.3.4 ECM vs. Cont at T28

Pa_28 <- extract_results(dds,vst,c(0,1,0,0,0,0,0,0,0,1),
                         default_prefix="Potri_ECM-vs-Cont_T28_",
                         labels=paste0(colData(dds)$Experiment,
                                       colData(dds)$Time),
                         sample_sel=colData(dds)$Time==28)
## There are 355 genes that are DE

4.3.5 Venn Diagram

4.3.5.1 All DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$all,
                            T7=Pa_7$all,
                            T14=Pa_14$all,
                            T21=Pa_21$all,
                            T28=Pa_28$all),
                       NULL,
                       fill=pal[1:5]))

4.3.5.2 UP DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$up,
                            T7=Pa_7$up,
                            T14=Pa_14$up,
                            T21=Pa_21$up,
                            T28=Pa_28$up),
                       NULL,
                       fill=pal[1:5]))

4.3.5.3 DOWN DE genes

grid.newpage()
grid.draw(venn.diagram(list(T3=Pa_3$dn,
                            T7=Pa_7$dn,
                            T14=Pa_14$dn,
                            T21=Pa_21$dn,
                            T28=Pa_28$dn),
                       NULL,
                       fill=pal[1:5]))

5 Session Info

## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] jsonlite_1.6.1              LSD_4.0-0                  
##  [3] limma_3.42.2                VennDiagram_1.6.20         
##  [5] futile.logger_1.4.3         forcats_0.4.0              
##  [7] stringr_1.4.0               dplyr_0.8.4                
##  [9] purrr_0.3.3                 readr_1.3.1                
## [11] tidyr_1.0.2                 tibble_2.1.3               
## [13] tidyverse_1.3.0             RColorBrewer_1.1-2         
## [15] hyperSpec_0.99-20200213     xml2_1.2.2                 
## [17] ggplot2_3.2.1               lattice_0.20-38            
## [19] here_0.1                    gplots_3.0.1.2             
## [21] DESeq2_1.26.0               SummarizedExperiment_1.16.1
## [23] DelayedArray_0.12.2         BiocParallel_1.20.1        
## [25] matrixStats_0.55.0          Biobase_2.46.0             
## [27] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
## [29] IRanges_2.20.2              S4Vectors_0.24.3           
## [31] BiocGenerics_0.32.0         data.table_1.12.8          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       rprojroot_1.3-2        htmlTable_1.13.3      
##  [4] XVector_0.26.0         fs_1.3.1               base64enc_0.1-3       
##  [7] rstudioapi_0.11        farver_2.0.3           bit64_0.9-7           
## [10] fansi_0.4.1            AnnotationDbi_1.48.0   lubridate_1.7.4       
## [13] splines_3.6.2          geneplotter_1.64.0     knitr_1.28            
## [16] Formula_1.2-3          broom_0.5.4            annotate_1.64.0       
## [19] cluster_2.1.0          dbplyr_1.4.2           png_0.1-7             
## [22] compiler_3.6.2         httr_1.4.1             backports_1.1.5       
## [25] assertthat_0.2.1       Matrix_1.2-18          lazyeval_0.2.2        
## [28] cli_2.0.1              formatR_1.7            acepack_1.4.1         
## [31] htmltools_0.4.0        tools_3.6.2            gtable_0.3.0          
## [34] glue_1.3.1             GenomeInfoDbData_1.2.2 Rcpp_1.0.3            
## [37] cellranger_1.1.0       vctrs_0.2.2            gdata_2.18.0          
## [40] nlme_3.1-144           xfun_0.12              rvest_0.3.5           
## [43] testthat_2.3.1         lifecycle_0.1.0        gtools_3.8.1          
## [46] XML_3.99-0.3           zlibbioc_1.32.0        scales_1.1.0          
## [49] hms_0.5.3              lambda.r_1.2.4         curl_4.3              
## [52] yaml_2.2.1             memoise_1.1.0          gridExtra_2.3         
## [55] rpart_4.1-15           latticeExtra_0.6-29    stringi_1.4.5         
## [58] RSQLite_2.2.0          highr_0.8              genefilter_1.68.0     
## [61] checkmate_2.0.0        caTools_1.18.0         rlang_0.4.4           
## [64] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [67] labeling_0.3           htmlwidgets_1.5.1      bit_1.1-15.2          
## [70] tidyselect_1.0.0       magrittr_1.5           R6_2.4.1              
## [73] generics_0.0.2         Hmisc_4.3-1            DBI_1.1.0             
## [76] pillar_1.4.3           haven_2.2.0            foreign_0.8-75        
## [79] withr_2.1.2            survival_3.1-8         RCurl_1.98-1.1        
## [82] nnet_7.3-12            modelr_0.1.5           crayon_1.3.4          
## [85] futile.options_1.0.1   KernSmooth_2.23-16     rmarkdown_2.1         
## [88] jpeg_0.1-8.1           locfit_1.5-9.1         readxl_1.3.1          
## [91] blob_1.2.1             reprex_0.3.0           digest_0.6.24         
## [94] xtable_1.8-4           munsell_0.5.0